Synchronization of globally coupled two-state stochastic oscillators with a 

state-dependent refractory period 
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We present a model of identical coupled two-state stochastic units each of which in isolation is 
governed by a fixed refractory period. The nonlinear coupling between units directly affects the 
refractory period, which now depends on the global state of the system and can therefore itself 
become time dependent. At weak coupling the array settles into a quiescent stationary state. 
Increasing coupling strength leads to a saddle node bifurcation, beyond which the quiescent state 
coexists with a stable limit cycle of nonlinear coherent oscillations. We explicitly determine the 
critical coupling constant for this transition. 



PACS numbers: 47.20.Ky, 47.54.-r, 82.40.Ck 



I. INTRODUCTION 

Many dynamical processes in nature can in the first in- 
stance be modeled as on-off processes. Examples include 
firing of neurons [1, 2], up-down spin systems, and blink- 
ing phenomena [3] . Although the dynamics of such on-off 
single systems is simple, when coupled together in large 
numbers they are known to exhibit complex global dy- 
namical behaviors such as synchronization, in which the 
entire system evolves as a single unit. The study of the 
synchronized dynamics of two-state networks is essential 
not only as a problem of fundamental scientific impor- 
tance [4], but also as a tool to help us understand many 
physiological processes [5]. For example, our cardiac sys- 
tem is a giant network of synchronized oscillators. A vast 
number of systems in physics and engineering also share 
such synchronization features [6, 7]. 

Coupled two-state units offer the simplest example of 
coupled units of discrete states. Such discrete models are 
useful for the study of phase synchronization of stochastic 
coupled oscillators [8, 9] and of excitable media [10, 11]. 
They can be invoked in ecological applications [12] where 
populations of various species interlinked through the 
food chain exhibit spatial and temporal synchronization. 
For sufficiently weak coupling between oscillators, the 
system as a whole typically reaches a stationary state 
where on average the dynamics ceases to exist. However, 
for relatively strong coupling each unit is driven away 
from its equilibrium, and the system as a whole may un- 
dergo a transition to a non-stationary state. 

Phase synchronization on networks has most com- 
monly been analyzed in terms of a Hopf bifurcation 
[8-11], that is, as a linear instability of the stationary 
state. Synchronization may also occur via other bifur- 
cation mechanisms such as, for instance, saddle node or 
homoclinic bifurcations [13]. A sub-critical Hopf bifurca- 
tion may, for example, involve a prior saddle node which 



can not be captured within the linear analysis. Such sub- 
critical scenarios have been observed in a variety of dif- 
ferent contexts including the Kuramoto model [14], Ku- 
ramoto model with time-delayed interactions [15], glob- 
ally coupled noisy bistable elements also with time de- 
layed interactions [16], and three-state coupled oscilla- 
tors [9]. These cases all involve the coexistence of a 
synchronous oscillatory state and a quiescent stationary 
state, with the corresponding hysteresis loop. 

The novelty of our model lies in the occurence of syn- 
chrony without destabilization of the stationary state, so 
that the oscillatory state may coexist with it without ex- 
hibiting hysteresis. In this scenario transitions between 
the synchronous phase and the asynchronous phase can 
only be induced by a strong perturbation. 

For a purely Markovian system the master equation is 
an ordinary differential equation. For coupled two-state 
systems it is a first order equation, which only has a fixed 
point. Therefore, in order to observe global synchroniza- 
tion in purely Markovian networks it is necessary to go 
beyond a two-state model [8, 9]. Converseley, to observe 
synchronization in a two-state system it is necessary to 
go beyond a simple Markovian model of coupled units. 
The synchronization of two-state units was perhaps first 
observed in the context of collective molecular motors, 
where spatial elastic coupling to a common environment 
leads to the occurrence of oscillations [17]. In globally 
coupled networks of two-state units, memory effects can 
also lead to synchronization [11]. In recent studies of cou- 
pled assemblies of genetic oscillators [18], memory (delay) 
effects have been found to play a significant role in the 
description of the observed synchronous dynamics of the 
on and the off (active and inactive) states of the system. 

Time delays in network dynamics can be introduced 
in many different ways. One is to introduce time de- 
lays in the interaction between the elements of the net- 
work [15, 16]. Another is to include a refractory period in 
the internal dynamics of the units [10]. Both of these to- 



gcthcr have also been considered [11]. Most of the models 
with a time delay that have been reported in the liter- 
ature consider uniform fixed time delays, be it in the 
internal dynamics and/or in the interactions. Huber and 
Tsimring suggested that temporally nonuniform (that is, 
distributed) time delays do not qualitatively affect the 
synchronization features observed in their model [16]. A 
more detailed analysis of the role of nonuniform time de- 
lays in synchronization has been reported in a cellular 
automaton version of the susceptible-infected-recovered- 
susceptible model [19], where a probabilistic refractory 
period is studied. 

In this work we also consider the effects of delay times, 
but of an altogether different kind than those described 
above. We consider identical two-state stochastic ele- 
ments with a fixed refractory period. The units are then 
coupled together all to all. The interactions among the 
elements affect the internal dynamics in that the refrac- 
tory period becomes dependent on the time-dependent 
state of the entire network. By combining numerical and 
analytical results, we find that a global dynamics may 
emerge as the coupling strength is increased beyond a 
critical value. We observe that synchronization occurs 
for some initial conditions but not for others. Unlike the 
usual case of a Hopf bifurcation, the synchronized behav- 
ior here emerges as a result of a saddle node bifurcation 
in the master equation phase space. 

The paper is organized as follows: In Sec. II we intro- 
duce our general model and show that no Hopf bifurca- 
tion is possible for this type of dynamics. In Sec. Ill we 
choose a specific model and show via numerical results 
that a global synchronized dynamics exists for our model. 
We also check the robustness of our results by considering 
some variants of the model. By computing a Poincare- 
type of map for the oscillatory state, in Sec. IV we de- 
mostrate that synchrony appears via a saddle node bi- 
furcation, and compute the critical point using this map. 
We confirm the critical point by direct integration of the 
generalized master equation. We end with some conclud- 
ing remarks in Sec. V. 



II. THE GENERAL MODEL 
A. The basic unit 

Our starting point is a binary unit characterized by 
two states, 1 and 2. Transitions from state 1 to state 
2 are governed by a Markovian (memoryless) process 
of rate 7, while transitions from 2 back to 1 occur af- 
ter a fixed delay time r measured from the time of the 
transition to 2 (Fig. 1). The waiting time distribution 
ipi^2{t) for a transition from 1 to 2 is thus an exponen- 
tial, = 7c -7 *, while that of a transition from 2 to 
1 is a (5-function, ip2^i = S(t — t). 

We next construct the generalized master equation for 
the probabilities Pi(t) to be in state i = 1,2 at time t. 
Normalization reduces this to a single probability, since 
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FIG. 1. Isolated binary unit. Transitions from state 1 to state 
2 occur at random times 7 -1 according to a Markov process 
of rate 7, while transitions from 2 to 1 occur at exactly time 
r after arrival at 2. 



Pi(t) = 1 - P 2 (t). Further, it is evident that 



(1) 



where J\^ 2 (t') dt is the probability of jumping from state 
1 to state 2 at some time within the interval [t\ t' + dt'}, 
with t' lying in the range \t — t, t] . (If t' were to lie be- 
low this range, the jump from state 2 back to state 1 
would have occurred already before time t.) In turn, 
J\^2 (f) dt' = jPi(t')dt', the probability of being in 
state 1 times the jumping rate to state 2. It follows that 



P 2 (t)= 1 dt' [I - P 2 (t')] 



(2) 



This linear integral equation for P 2 can also be written 
in differential form to arrive at the master equation 



CP 2 = IT, 



(3) 



where the linear operator C involves infinite derivatives, 



(-rf 
(n + 1)! \dt 



(4) 



The master equation (3) predicts that, after a tran- 
sient, the system reaches a stationary equilibrium, 
Ps{t) P^i which is determined by the rate 7 and the 
refractory period r via the equation 



P2 



1 T 



1 + 7T 



(5) 



This result is easily understood: at equilibrium, the mean 



time spent by the unit in state 1 is t\ 



7" 



while 



the mean time spent in state 2 is exactly the refractory 
period, T2 = t. The equilibrium probability P 2 * is then 
simply the fraction of the time spent in state 2, 



P* 



T 2 



Tl +T 2 



(0) 



which immediately leads to Eq. (5). 

We can furthermore see, as follows, that the equilib- 
rium probability P 2 * is always the single stable solution. 
Equation (2) is linear and invariant under time transla- 
tion. We can therefore pick any time to and wc know that 



3 



N 2 (t)/N °- 7 
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FIG. 2. Ensemble of TV = 2 x 10 4 isolated units with 7 
and t = 1, which leads to P 2 * = 1/2. 



if a unit is in state 2, it must have transitioned there at 
some time during the time interval (to — T,to). Had it 
transitioned to state 2 earlier than that, say in the in- 
terval (to — 2t, to — t) , it would have returned to state 1 
before time to- If h transitioned to 2 before time to — 2r, 
it might return and be back in state 2 at time to, but 
this event is already accounted for in the accounting of 
transitions in the interval (to — t, to). One can thus write 
the transient solution as 



(7) 



where the set {qn}^ = i depends on the "initial condition," 
that is, on the function P2{t) in the time interval [to — 
r, to], and {An}^^ is the set of solutions of the equation 
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(8) 



The real and imaginary parts separately, A„ = p n + iuj n 
obey the equations 



e P " T cos (w n r) 
e~ PnT sin (u) n r) 



1 = 7 Pn 



(9) 



From Eq. (8) it is clear that the solution p n = uj n = 
is spurious. Furthermore, the choice p n > leads to an 
immediate contradiction, namely, that 7 < 0. This is 
impossible since 7 is a rate, hence positive. Therefore 
the real parts of the A„ must all be negative, and thus 
P 2 {t) goes to P 2 for any initial condition. Indeed, the 
dynamics of an ensemble of N such isolated units consists 
of a transient of damped oscillations and eventual arrival 
at equilibrium, with fluctuations due to the finite number 
of units. A typical evolution of N such independent units 
is shown in Fig. 2, where we show P2(i) = N2(t)/N vs 
t. Here N 2 (t) is the number of units in state 2. The 
small fluctuations just visible in the figure can be further 
mitigated by working with a larger number of units. 



B. Coupled array: Generalized master equation 

Next we globally couple N of these two-state oscilla- 
tors, that is, each oscillator is coupled to all the other 



oscillators. We point with special emphasis to the nature 
of the coupling about to be described, because it is this 
feature that distinguishes our model from others in the 
literature. 

We take the transition rate from state 1 to state 2, and 
the refractory period in state 2 to become functions of the 
global state of the system. These are therefore themselves 
now time dependent. That is, if at a time t there are N\ 
units in state 1 and N 2 units in state 2 (always with the 
constraint N = N\ + N 2 ), the time-dependent transition 
rate and refractory period are modified from their values 
in the isolated units, 



7 -> lit) = 7 
r — > r(t) = t 



N 2 (t) 

N 
N 2 {t) 

N 



70P 2 (t)), 
r(P 2 (i)). 



(10a) 
(10b) 



Here we have indicated the implementation of the ther- 
modynamic limit N — > 00, 



N 2 (t) 
N 



(ii) 



and we recall again that Pi(t) = 1 — P 2 (t). We have not 
yet specified the explicit form of the dependences on the 
state, but only that the transition rate from 1 to 2 and 
the transition delay from 2 back to 1 are now functions 
of the state. Recall that in the single units, an oscillator 
that arrived at state 2 at time t — r jumped back to state 
1 at time t. In an uncoupled collection of single units, 
the times at which different units are observed to leave 
state 2 are random because the arrival times are random. 
However, the length of time spent in state 2 after each ar- 
rival is fixed at the value r. In a coupled array once again 
units arrive at state 2 at random times but remain there 
for a fixed time period. Because of the random arrival, 
the departures are also randomly distributed. However, 
there are now two differences. The rate at which units 
leave state 1 is not fixed but depends on the instanta- 
neous state of the coupled system. Similarly, the waiting 
time in state 2 varies with system state. We clarify what 
we mean by the time dependence of the refractory period, 
without yet fixing the form of the dependence implicit in 
Eq. (10b): all the oscillators that arrived at state 2 at 
time t — r(t) jump back to state 1 at time t. Thus if 
one has a number of units in state 2 at time t, some will 
jump back to state 1 right then and some not, depending 
on their arrival time. We emphasize that this is not a 
stochastic feature of r(t) but rather a consequence of the 
stochasticity in the arrival times embodied in 7(f). 

To accommodate these time dependences in a mean 
field description appropriate for all-to-all coupling, we 
generalize the master equation (2), 



W) = f 

J — ; 



dt'F (P 2 (t'))Q(t,t'). 



(12) 



Here 



F (P 2 (t')) dt' - 7 (iMO) [1 - p 2(t')] dt' (13) 
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is the probability to jump from state 1 to state 2 at some 
time within the interval [t' , t'+dt']; the factor 1 — P 2 (t') = 
Pi(t') is the probability to be in state 1 at time if, and 
7 {Piit')) is the (ever-changing) rate at which the jump 
from 1 to 2 takes place at time t' . The factor 0(t, t') keeps 
track of this arrival rate and "freezes" each oscillator in 
state 2 for a precise length of time (refractory period) 
determined by the state of the system when the oscillator 
arrived there. 

We have highlighted this last term because the core 
and novelty of our model lie in it. Consider the units 
that arrive at state 2 at time t' . The number of arrivals 
is stochastic since the departure from state 1 is random. 
Furthermore, the rate of random arrivals depends on the 
state of the system at that time, 7 (P 2 (<')). The depar- 
ture of these units from state 2 occurs at a definite delay 
time t" later, t" —t' = r(t"), which depends on the state 
P 2 {t") at the time of departure. If t" - t' < r(<") we 
know that the units that arrived at if have not departed. 
On the other hand, if t" —t'> r(t") we know that they 
have. Thus, if the delay time r(t") is not reached at any 
time t" in the interval [if, t], we have that 8 (t, t') — 1, an 
indication that all the oscillators that arrived in state 2 
at time t' are still there (while others may continue to ar- 
rive). If the delay time is reached at any time within the 
interval, then Q(t,t r ) = because the units that arrived 
in state 2 at time if have left again. This step-like func- 
tion character arises from the fact that the delay time is 
fixed by the population probabilities, that is, it is not a 
random variable. Explicitly, © is the step-like function 



Specifically, such a stationary solution would satisfy 



e(M') 



1 ifVt" G [t',t], t"-t' <r[P 2 {t")} 
otherwise 

(14) 

We call the condition that leads to Q(t,t') = 1 the 
inequality. In Appendix A we explore this question in 
detail independently of the particular functional form of 
7 and t as functions of time. Here we take advantage 
of the lessons learned there and implement the behaviors 
that arise in our case. 

Our ultimate goal is to find the behavior of the Pi(t) 
at long times. In particular, we arc interested in estab- 
lishing whether this behavior is time dependent or time 
independent, and under what conditions one might ob- 
serve the former or the latter. 



C. Stationary solutions 

We start by exploring the possible occurrence of one 
or more stationary states, that is, one or more time- 
independent solutions P 2 (i) = P 2 *. This case is covered 
by the first part of Appendix A and falls in the class 
of behaviors described by the generalized master equa- 
tion (A2), 



Jt-T 



(P2(t)) 



dt'F{P 2 {t')) 



(15) 



P* = 



l + 7 (P 2 *)r(P 2 *)' 



(16) 



Note that this is only similar to Eq. (5) for the uncoupled 
system in aspect since Eq. (16) is in general a non- linear 
equation for quiescent states. This nonlinear equation 
could have more than one solution (multi-stability) or 
even no solution at all, depending on the particular func- 
tional forms chosen for 7 (P 2 ) and r (P 2 ). 



D. Absence of Hopf-type instabilities 

While we have not yet established whether Eq. (16) 
actually has one or more solutions, we can show that 
the coupled array does not support Hopf-type instabili- 
ties leading to oscillatory or chaotic solutions. Indeed, a 
linear analysis at best leads to other quiescent solutions. 
To show this we introduce a perturbation, 



P2it) = P 2 +ee 



j.l I 



(17) 



with |e| < 1. First, we note that if r'(P 2 (i)) ~ 0(1), as 
it will be in our model explicitly introduced later, then 
dr (P 2 {t)) jdt ~ e < 1. (The prime denotes a deriva- 
tive with respect to the argument.) According to the 
discussion in the appendix, the generalized master equa- 
tion (15) is valid if fi is negative and for sufficiently short 
times even if \x is positive. To first order in e, the master 
equation is then given by 

ee"* =e / dt'F' (P 2 *) e' 1 *' 

+ / dt'F(P*). 

Jt-r(p;)-er'(PS)e^ 

Straight substitution of Eq. (16) in this equation shows 
that /i must satisfy 



e-» T - 1 



(18) 



where 



t*=t(P 2 *) 



r = 



P'(P*) 



t> (P 2 *)F(P 2 *)-1 



Therefore, all the information about the linear stability 
of P 2 * resides in the two numbers t* and T. No additional 
information about the specific functional form of t(P 2 ) 
or of 7(P 2 ) is required. 

To analyze equation (18) one can separate the real and 
imaginary parts, (j, = J3 + iQ, in terms of which condi- 
tion (18) can be separated into the two equations 



1 



cos (fir*) = T" 1 /?, 
sin (fir*) = r^ft. 
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A Hopf-typc instability requires that a stable stationary 
solution become unstable with non-zero frequency. A 
critical point must satisfy (3 C = with Q c ^ 0, that is, 

cos(f2 c T*) = 1 and sin(f2 c T*) = r _1 f2 c . 

However, this pair of equations only has the solution 
il c = 0. Note that the solution [i = is no longer nec- 
essarily spurious as it was in the decoupled array. There 
the solution A„ = required that ■yr = — 1 [cf. Eq. (8)], 
which is impossible because both parameters are positive. 
However, here the condition IY* = — 1 can be seen as an 
equation for a critical value of the control parameter. 

In summary this linear analysis predicts no complex- 
ity beyond possible bifurcations between stationary so- 
lutions. However, the analysis itself is limited: our non- 
linear infinite-dimensional system has the potential to 
exhibit more complex behavior such as oscillatory or 
chaotic orbits. We will next explore the possible occur- 
rence of limit cycles (synchronization of the units) via 
global bifurcations, that is, bifurcations that go beyond 
the vicinity of a fixed point to involve larger regions of 
phase space. For low dimensional systems there are three 
generic global bifurcation pathways to cycles [13] : saddle- 
node, infinite-period and homoclinic bifurcations. We 
find below that our infitine-dimcnsional system also fol- 
lows one of these pathways. 

E. Generic features beyond stationary states 

We continue to explore possible behaviors, now away 
from stationary states. Our analysis in Appendix A 
shows that there is only a limited array of possible be- 
haviors. We will basically encounter only two. 

The first, which we will call behavior B, occurs when 
S(t,t') can be expressed as a Heaviside function, cf. 
Eq. (Al) in the appendix, 

Q(t 1 t') = 6(t' -t + r(P 2 (t))). (19) 

As shown in the appendix, this expression is valid when 
(1) r(t) is continuous and piccewisc differentiable; (2) 
where differentiable, dr/dt < 1; (3) r(t) may have discon- 
tinuities presenting a downward jump from a high value 
to a low value as time increases, but not from a low value 
to a high value. 

The other behavior we will encounter, behavior A, oc- 
curs in one of two ways. It can come about because (3) 
above is violated, that is, r(t) exhibits a discontinuity 
presenting an upward jump from a low value to a high 
value, or because condition (2) above is violated. There 
are several ways to violate the condition dr/dt < 1 and 
thus to enter regime A. These four ways are illustrated 
in Fig. 3. One, illustrated in the upper right panel of 
the figure, occurs smoothly. It starts at time t\ where 
dr/dt\ t=tl = 1. A second way to enter regime A is il- 
lustrated in the lower left panel. Here regime A again 
begins at time t\, but r(t) need not be differentiable at 




t, h t, h 



FIG. 3. Four ways to enter regime A from regime B. Upper 
left: when there is a discontinuity in r from a low value to a 
high value. Upper right: when dr/dt smoothly goes through 
the value dr/dt = 1 from whence it increases. Lower left: 
similar to upper right except that r(t) need not be differen- 
tiable at the point of entry in A. Lower right: here there is a 
discontinuity in r from a low value to a high value together 
with a derivative that approaching the point of discontinu- 
ity that exceeds unity, lim + (dr/dt) > 1. In all four cases 
entry into regime A occurs at time t\ and exit at time ti. 

t\. Instead, here lim t t + (dr/dt) > 1. The third way to 
enter regime A, illustrated in the lower right panel, is via 
a discontinuity from a high value to a low value at time 
t\, but with a violation of condition (2) above, that is, 
again lim. .+ (dr/dt) > 1. The final way, illustrated in 
the upper left panel, occurs when there is a discontinuity 
from a low value to a high value. In all four cases to end 
regime A and re-enter behavior B it is not sufficient for 
dr/dt < 1 for some t > t\. The end point occurs when 
r(t) crosses the line t — ti + r(ti), that is, at time t 2 that 
satisfies 

T(h)=t 2 -t 1 + r{t 1 ). (20) 
The associated O then is 

e(t,t') = 9(t' -h+rih)), tG[ti,t 2 ]. (21) 

We conclude with a comment: we have indicated that 
in all cases regime A ends with a re-entry into Regime B 
in a finite time. This is, in fact, not necessary. Regime A 
could continue forever. But this would imply a continual 
growth of r(i), which leads to the uninteresting conclu- 
sion that (since exit from state 2 of our array units is 
forever slowed down) all units will end up in state 2 and 
will remain there. This is a trivial dynamics that will not 
be pursued further. 

We end this section by noting that we have so far 
made no specific assumptions about the functional forms 
of r y(P 2 (t)) or r(P 2 (t)) but have, instead, discussed the 
problem in full generality. We also note that our discus- 
sion so far is fully analytic. Now, to proceed we introduce 
specific forms for these functional dependences. 
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III. THE SPECIFIC MODEL: 
SYNCHRONIZATION 

A. Choosing a model 

At this point any further analysis of the emergence of 
synchronization in our ensemble of globally coupled oscil- 
lators requires us to specify the time-dependent quanti- 
ties 7(i) and r(i). For the rate 7 wc choose the standard 
(Arrhenius-type) exponential law [8-11] 



7(t) = gexp 
= gexp 



a[N 2 (t) - mjt)} 
N 

a[2N 2 (t) - N] 
N 



(22) 



For the refractory period we introduce a non-monotonic 
dependence on N2(t)/N, 



"(*) = ^0 



N^N^t) N 2 (t)(N - N 2 (t)) 



N 2 



tq- 



N 2 



(23) 



which exhibits a maximum when half the units are in 
state 1 and the other half in state 2. An imbalance in 
either direction shortens the refractory period. Here g, 
a and tq are coupling parameters. Next, we can scale 
time so that with no loss of generality one of the param- 
eters can be fixed, say g = 1, and hence the system has 
a total of three free parameters: a, tq and N. We fo- 
cus on the thermodynamic behavior N ^> 1, so that the 
phase space diagram of the system is, in fact, the two- 
dimensional parameter plane (a, To). We have observed 
that varying just one of these two parameters allows us 
to capture the main features of the system (the transi- 
tion to synchronization), so we fix To = 2 and present our 
analysis as a function of the coupling parameter a. 

Note that in principle the coupling parameter may be 
positive or negative. A positive parameter implies that if 
many units are in state 1, say, then they leave that state 
more slowly than if there are only a few. We choose to 
study the case of negative a, where the behavior is the 
opposite, namely that units leave state 1 at a greater rate 
when it is highly occupied. That is, the system "avoids" 
crowding. This is consistent with the form chosen for the 
refractory period: it, is shortened when either state is 
overcrowded. 



B. Outcomes in the thermodynamic limit 

First we test the mean field theory for our model, that 
is, we solve Eq. (12), which of course must be done nu- 
merically. We observe that after a transient dynamics, 
the system approaches a long-time behavior that depends 
on the choice of the coupling parameter a. When \a\ 
is small ("small" to be specified accurately later), the 
system approaches a quiescent state well described by 
Eq. (16). On the other hand, when the coupling is suffi- 
ciently strong (to be determined later) , the system either 



approaches this quiescent state, or, following a transient, 
the system settles into a globally synchronized dynamics 
wherein the probabilities Pi(t) oscillate in highly anhar- 
monic but essentially periodic fashion. Fig. 4(a) displays 
a typical such global limit cycle. We thus observe the 
coexistence of a stationary state and a stable limit cy- 
cle. Note that wc can not necessarily simply search for 
initial conditions P 2 {0) that lead to one behavior or the 
other. Because the system has a memory, the "basins of 
attraction" in general would have to be characterized in 
functional space and require knowledge of P 2 (t) for past 
times before t = 0. For our model, we need to know P 2 (t) 
at most in the interval [— to/4, 0]. If all the units are ini- 
tially in state 1, that is, if ^(0) = 0, then we do not need 
information about the past. We have not attempted to 
unravel the functional initial condition problem here. 

We test the mean field master equation by also di- 
rectly simulating the equations of motion of our system. 
Our simulations here are carried out for N = 10 4 . Here 
again we find that for sufficiently strong coupling, after 
a transient dynamics the system approaches one of two 
long-time behaviors: the system either converges to a 
stationary state wherein Ni and N 2 approach essentially 
constant values (again well described by the mean field 
quiescent state (16)), or the system goes to a final state 
of synchronized dynamics that consists of coherent oscil- 
lations of Ni and N 2 . The number N = 10 4 is sufficiently 
large for us not to observe fluctuations that would per- 
turb the system out of an oscillatory dynamical behavior 
during the time of our simulations, but presumably a suf- 
ficiently long run would exhibit such transitions. In any 
case, we thus again observe the coexistence of a station- 
ary state and a stable limit cycle. Fig. 4(b) displays a 
typical such global limit cycle, this one for the same ini- 
tial condition and coupling as in Fig. 4(a). We see that 
N 2 (t) /N exhibits very anharmonic oscillations with well 
defined amplitude and frequency. The similarity between 
the two panels reassures us that N = 10 4 is sufficiently 
large to capture all the essential features of the thermody- 
namic behavior of the system. The excellent agreement 
between the solution of the mean field master equation 
and the direct numerical simulations is gratifying. The 
simulations inevitably exhibit very small fluctuations not 
present in the mean field model because the simulations 
are carried out for finite, albeit large, N. 

We will subsequently argue that the two solutions, syn- 
chronous and asynchronous, are co-existing stable solu- 
tions, and that the synchronization is here not related 
to the destabilization of the stationary state. Regardless 
of initial condition, a sufficiently strong external distur- 
bance applied to the system can drive it from one state 
to the other. In a system with a finite number of units 
the two solutions are metastable states not impervious to 
fluctuations; indeed, a sufficiently strong fluctuation can 
cause a transition from one of these states to the other. 
Thus, now the perturbation need not be external to the 
system. 

In the next subsection we present a number of results 
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P* ~ 0.422. 
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FIG. 4. (Color online) Long-time behavior of the model for 
a = —2 and to = 2. (a) Numerical solution of the mean 
field equation (12). (b) Direct numerical simulations of the 
equations of motion for TV = 10 4 units coupled according to 
the prescriptions (22) and (23). 

that arise from direct numerical simulations of coupled 
units. Then, in the next section, we return to the master 
equation and, to understand the synchronization proper- 
ties from this more analytic point of view, we construct 
a Poincarc map for the oscillatory state. Among other 
results, we calculate the critical coupling for our model 
in the thermodynamic limit. 

C. A variety of results of direct numerical 
simulations 

We present several sets of results, not all meant to be 
quantitatively accurate but rather to illustrate a num- 
ber of interesting observed features of our system. More 
detailed studies of these features will be explored in sub- 
sequent work [21]. 

First, we illustrate the fact that for small N the fluc- 
tuations are sufficiently strong to drive the system back 
and forth between oscillatory and stationary behavior. 
We have not explored the dependence of the rate of these 
transitions on N, but note that they arise from the mi- 
croscopic dynamics. Three realizations of such a progres- 
sion in time are shown in Fig. 5 for an array of iV = 100 
all-to-all coupled units. The precise parameter values are 
not particularly important other than to confirm that the 
stationary episodes in the figure arc consistent with the 
stationary value obtained from Eq. (16), which here is 
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FIG. 5. (Color online) Typical realizations for an array of 
100 units coupled according to the prescriptions (22) and (23) 
above critical coupling. The parameter values are the same 
for the three realizations. The mean field stationary value for 
these parameters is P 2 * — 0.422. 

A second point to explore is the robustness of the re- 
sults to variations in the model. Ours would not be a 
very interesting effort if it depended strongly on the spe- 
cific features of the model without room for variation. 
We have explored two specific variations. 

To display the first generalization we note that we can 
express our basic model more generally than we have 



done, namely, we can introduce a rate gi(t,t') of transi- 
tions from state 1 to state 2, and likewise a rate g 2 (t,t') 
for transitions from state 2 back to state 1. Here t' is 
absolute time, and t denotes the waiting time. We con- 
tinue to assume that transitions from 1 to 2 constitute 
a Markov process so that gi(t,t') = j(t') as indicated in 
Eq. (22). However, for the return we allow a distributed 
waiting time of the form 



92 t, 



1 



t 



N J T \T{N 2 (t')/N) 



(24) 



which reduces to our working model with a single waiting 
time when b — > oo. The full master equation is 



P 2 {t) = / rfi" 7 (P 2 (t - t")) (1 - P 2 (t - t")) 
Jo 



x cxp 



df 



f 



o to \r(P 2 (t-t" + t')) i 



(25) 



which reduces to Eq. (12) with (14) when b — > oo. We 
have ascertained for various values of b that this also 
leads to oscillatory solutions for sufficiently strong cou- 
pling, that these coexist with stationary solutions in this 
regime, and that the solutions obtained from the sim- 
ulations of a large number of coupled units and those 
obtained from the generalized master equation coincide. 
The stationary solutions of the master equation satisfy 



C 7 (P 2 *)[t 2 (P 2 *)]^ 
1 + C 7 (P 2 *)[t 2 (P 2 *)P 



(26) 



with 



C= [r (b+l)]*Fir 



6 + 2 

b + i 



Here T is the gamma function. It is straightforward to 
verify that this reduces to Eq. (16) when b — > oo. The os- 
cillatory solutions exhibit the same strongly anharmonic 
behavior for the values of b we have explored (in the 
neighborhood of b — 100) as they do in our working 
model. A typical realization of such a stationary solu- 
tion is shown in the upper panel of Fig. 6. 

We have also noted that P 2 = is a stationary solu- 
tion of our working model as well as of the generalization 
discussed above. To vary this, we have explored the con- 
sequence of modifying the waiting time in state 2 from 
the form r(t) of Eq. (23) to r(i) + Stq. We have tested 
the effect of this addition for small values, Stq ~ 0.1, and 
have again found the same behaviors. A typical realiza- 
tion of an oscillatory solution for this model is seen in 
the lower panel of Fig. 6. 

Finally, we illustrate in Fig. 7 the loss of synchroniza- 
tion with our working model as we cross the critical cou- 
pling. The simulation in the figure is for 160, 000 coupled 
units. In panel (a) the coupling is well above the critical 




FIG. 6. (Color online) (a) Typical realization for direct sim- 
ulations of the model described in Eqs. (24) and (25) with 
ro = 2, a = —4, b = 100, and N = 10, 000 coupled according 
to the prescriptions (22) and (23). (b) Typical realization for 
model with "displaced waiting time" with parameters to = 2, 
a = -5, 5t = 0.1, and N = 10, 000. 



value; in (b) it is just below the critical value, and in (c) 
it is well below this value. 

As promised, we now return to the master equation 
and, to understand the synchronization properties from 
this more analytic point of view, we construct a Poincare 
map for the oscillatory state so as to obtain a quasi- 
analytic estimate of the critical coupling in the thermo- 
dynamic limit. 



IV. POINCARE- TYPE MAP FOR THE 
OSCILLATORY STATE 

To understand the synchronization properties of our 
coupled array, and guided by our numerical solutions but 
proceeding as analytically as possible, we focus on a sin- 
gle typical oscillatory cycle, cf. Fig. 8(a). We define a 
cycle as the evolution between the two evident disconti- 
nuities in the figure, indicated by dashed lines. In more 
detail, we observe that at the start of a cycle, P 2 (t) at first 
grows gradually and then it quickly decreases to a mini- 
mum value that is in general different from the starting 
value of the cycle. Then the process begins again. The 
cycle shown in Fig. 8 is a generic n-th cycle, where t n 
denotes an arbitrary time which does not play any ex- 
plicit role in the calculations. The n-th cycle starts from 
some initial condition P 2 (t n ) = p n which is the end point 
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FIG. 7. Numerical simulation of N = 160, 000 units: (a) well 
above critical coupling; (b) just below critical coupling; (c) 
well below critical coupling. 



of the previous discontinuity (and whose value depends 
on the past). After a time T 2 the system undergoes an- 
other discontinuity ending at time i„+i = t n +T 2 , that is, 
P2{t n +i) = Pn+l- Note that T\ and T 2 change from cycle 
to cycle because periodicity is not perfect, that is, they 
both depend on n via p n . The key to the synchronization 
is the Poincare-type map 



Pn+l = / (Vn) ■ 



(27) 



The oscillatory solution corresponds to a fixed point of 
the map, p* = f(p*), and it is stable if \f (p*)\ < 1. 
Otherwise it is unstable [13]. 

In order to solve the master equation (12) and com- 
pute the map (27), we note that each cycle exhibits two 
significantly different regimes of behavior, indicated as A 
and B in Fig. 8. These are precisely the regimes A and 
B discussed in detail in Appendix A and in Sec. II. 



FIG. 8. (Color online) (a) P2-profile during the n-th cycle, 
(b) Transition delay time r vs time (solid, blue), and the 
straight line r = t — t n + t(p„) vs time (dashed, red). The 
A portions use the analytic result Eq. (31), while the results 
for regime B are outcomes of the numerical integration of 
Eq. (34). Parameter values: to = 2, a — —1.5 and p n — 0.05 
(pn+i was inserted manually in order to clarify the plots). 



A-regime: This regime is characterized by the fact 
that the curve T(P 2 (t)) lies above the line t — t n + r(p n ), 
t (P 2 (t)) >t-t n + r(p n ), for t n <t< t n + Ti. 

B-regime: Here the inequality is reversed (the curve 
lies below the line), r (P2W) < t-t n +r(p n ) 1 for t n +Ti < 
t<t n + T 2 . 

We next consider the generalized master equation in 
each regime separately. Consider first regime A. From 
Eq. (21) it follows that for t within the interval [t n ,t n + 
Ti], 



e(t,t') = e(t' -t n + T(t n )), 



(28) 



where for economy of notation we have written 
r(P 2 (in)) = T (t n )- The master equation can then be 
written as 



P2 A (t) 



dt'F(P 2 (t')) 



(29) 



lt n -r(t n ) 

We can alternatively write this in the differential form 

(30) 

with the initial condition P 2 (t n \p n ) = p n . The solution 



P 2 A = F(P 2 A ) 
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(*;Pn) = 1 - T^-ei" 1 (ei (2a (1 - p n )) - e a (t - t n )) 



2a 



where ci denotes the exponential integral, 



(31) 



ei(z) 



du- 



and ei -1 denotes its inverse. Hence, in this regime we 
have presented an exact analytic solution. 

We carry this analytic focus further by noting that the 
A- regime ends after a time interval Ti, which is defined 
by the equation 



Ti+t (p n ) = t (P 2 A (Ti + t n ;p n )) , 



(32) 



which identifies the point at which the r curve and the 
line t — t n + r(p n ) intersect. This identification is im- 
portant because reaches a maximum at this point 
(followed by a plateau discussed below), and with the an- 
alytic solution (31) we can compute this maximum and 
the time at which it occurs. Note again that Ti depends 
on p n , but does not depend on t n . 

In the B-regime we are forced to rely partly but min- 
imally on numerical simulations, which lead to the ob- 
servation that dr (P 5 (t;p n )) /dt < 1. Here we have ex- 
plicitly indicated the dependence of r on t via P 2 B . The 
master equation in this regime can be written as 



dt'F{P 2 {t'-p n )) (33) 



(cf. Eq. (A2)). Taking a derivative of this equation with 
respect to t and solving for P 2 (t]p n ) yields the differen- 
tial form 



where 



V(t) 



F(p 2 B (t lPn ))~F(r (t- T (p 2 B ))) 

l-F(V(t-r(P 2 B )))r'(Pf) 



P 2 B (i;p„_i) ifi<i n 
P 2 A {t;p n ) iit n <t<t n +T l 
P 2 B (t;p n ) iit>t n +T 1 



, (34) 



together with the initial condition 

P 2 S (Ti + t„;p n ) = P 2 A (Ti + t n ;p n ) 



(35) 



Equation (34) is a non-autonomous differential equation 
for t— t (P 2 B (t;p n )) < t n because it depends on P 2 in the 
previous cycle, P 2 (t\p n -\). It is also non-autonomous 
for t n < t — t (P 2 B (t;p n )) < t n + Ti because it de- 
pends on P 2 in the same cycle, P 2 {t;p n ). In general, 
for t — t (P 2 B (t;p n )) > t n + Ti it is a delay differen- 
tial equation [20] because it only P 2 at previous times; 
however, as seen immediately below, this regime is not 
reached in our system in the oscillatory state. 



More specifically for our case, the description of the 
system in regime B then proceeds as follows. For t — 
t (P 2 S (t;p n )) < t n , the P 2 -profile exhibits a plateau as 
seen in Fig. 8(a). Note that this is a word used here to 
describe a numerical observation rather than an analytic 
precise outcome. Then, for t n <t — r (P 2 (t;p n )) < t n + 
Ti it quickly decreases up to the time at which it exhibits 
the discontinuity. This occurs at t n + T 2 , which satisfies 
T 2 — r (P 2 B (t n + T 2 ;p n )) < Ti. Equation (34) remains a 
non-autonomous differential equation during this entire 
portion of the cycle, that is, it reaches the discontinuity 
before it becomes a delay differential equation. 

At the level of Eq. (34) the discontinuity is related 
to the divergence Pf (t — >• t n + T%\p n ) — > —00. At the 
discontinuity, P-zif) falls to the value p n +i, and regime A 
begins again, 



Pn+l — (t n+ i, p n+ i) . 



(36) 



We also note that numerical simulations show that p n is 
small, so that we can assume that T 2 — r(f„ + i) > T±. In 
this case, t lies in the range, [t n + T 2 — r(t„_|_i), t n + T 2 ]. 
This places us in regime B of cycle n. We can thus write 



Pn 



+ r(p„ + i) 



dt'F(P 2 B (t';p n )) (37) 



where t n+ \ = t n +T 2 and where we have invoked Eqs. (37) 
and (29). This is an equation that relates p n +i to p n , that 
is, it implicitly contains the map (27). 

To compute the map we have again combined nu- 
merical with analytical results by numerically solv- 
ing Eqs. (34) and (37) using the analytic solution for 
P 2 A {t\p n ). Our results for r = 2 are summarized in 
Fig. 9, which clearly predicts that synchronization ap- 
pears via a saddle node bifurcation. For a > a c the 
function / has no fixed points, so the system is not able 
to synchronize. At a = a c the function / just touches the 
line Pn+i = Pn, and for a < a c there appear two fixed 
points related to periodic orbits. The lower fixed point 
corresponds to the stable oscillation that we observe in 
our simulations; the other fixed point corresponds to an 
unstable limit cycle. From this analysis we find that the 
critical point for to = 2 occurs at a c = —1.42. 

We can also estimate the critical point directly from 
the integral equation (12). To arrive at an ever more 
accurate critical point via this route, we must use ever 
smaller time steps At in the numerical integration. In 
Fig. 10 we present the results as we decrease the time 
step, and obtain the value a c = —1.44 upon extrapolation 
to At — > 0. The two methods to calculate the critical 
point thus show a high level of agreement. 



V. CONCLUSIONS 

It is well known that a globally syncrhonized state can 
not occur in Markovian arrays of coupled two-state sys- 
tems. A number of models of coupled arrays of two-state 
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FIG. 9. (Color online) Computation of the map (27) for to 
2 and three different values of a. 




FIG. 10. (Color online) Variation of the critical point with 
the integration step size At (see text). The continuous curve 
is the best power law fit to the numerical data, —a c — 
0.303(At) ' 846 + 1.44. Extrapolating to At -> leads to 
the critical point a c — —1.44, which is close to the value 
a c — —1.42 obtained from the saddle node analysis of the 
master equation. 



nized dynamics, and that incorporates a physical fea- 
ture endemic to two-state systems: a refractory period. 
Specifically, in our model the single units are not Marko- 
vian to begin with but instead, upon arrival at one of the 
two states, they are subject to a delay time before they 
can return to the other state. The coupling among units 
is reflected in a state dependence of this delay time. In 
other words, the refractory period itself is affected by the 
state of the system and therefore changes with time. 

In this globally coupled array we have found a synchro- 
nized state that exhibits highly anharmonic dynamics 
with cyclical discontinuities. The global dynamics of the 
network can be described by a generalized master equa- 
tion. We have solved the generalized master equation to 
show that this type of synchronization is the result of a 
saddle node bifurcation in the phase space of the system. 
The critical coupling strength, as well as the amplitude 
and period of oscillations, are derived analytically by ex- 
plicitly solving the generalized master equation. We find 
that as we approach the critical coupling parameter a c 
from below, transient oscillatory behavior of longer and 
longer duration T occurs, with T —> oo as a — > a c . We 
are currently exploring the functional dependence of T 
on \a — a c \. 

The main points to be emphasized can be summarized 
as follows. Incorporating a refractory period that de- 
pends on the global state of the system has given rise 
to a rich dynamics, one that can exhibit oscillations and 
coexistence between the synchronous and asynchronous 
states (bistability). The nature of the synchronization 
that leads to self-organization far from equilibrium is here 
not related to the destabilization of the thermodynamic 
branch. Indeed, this branch, which leads to stationary 
behavior, remains stable and coexists with the stable os- 
cillatory branch. This is an entirely different synchro- 
nization paradigm than that obtained from the usual 
Hopf bifurcation and the attendant destabilization of the 
thermodynamic stationary state. We have ascertained 
that our model is robust by considering two variations 
of the model and showing that all three models lead to 
similar synchronization properties. 

Refractory properties are evident in a large variety of 
systems ranging from the physiological (refractory prop- 
erties of an excitable membrane) to the physical (blinking 
quantum dots). The rich dynamics of our simple on-off 
model makes this an excellent candidate for adaptation 
to specific applications. 



units have been constructed that depart in one way or 
another from Markovian dynamics and that as a result 
do lead to global synchronization. These models share a 
feature: the memory lies in the coupling between units 
and does not affect the dynamics of the units themselves. 

We have here constructed a novel departure from 
Markovian behavior that also leads to globally synchro- 
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Appendix A: Generalized master equation for mean 
field coupling 

In this Appendix we explore possible explicit func- 
tional forms of r and the concomitant behavior of 9. 
Note that we are ultimately especially interested in find- 
ing the long-time behavior of the coupled array. To make 
the case clear, we build up from the simplest possible sce- 
narios to the most complex that we will encounter. 

Consider first the case of a function r(P2(i)) that de- 
creases smoothly with time, as sketched in the upper left 
panel of Fig. 11. The construct in the sketch clearly 
shows that the inequality is valid only at the foot of the 
triangle. Units that arrived earlier than t — r(t) will have 
left before time t and thus do not obey the inequality. 
The function can then be expressed in the simple form 

e(t,t') = e(t' -t + T (p 2 (t))), (Ai) 

where 9 is the Hcaviside step function. The resultant 
master equation reads 

P 2 (t)= [ dt'F(P 2 (t')). (A2) 

Jt-T(P 2 {t)) 

If r(t) increases smoothly with time, then the same anal- 
ysis is valid provided dr/dt < 1. This is apparent in the 
schematic in the upper right panel of the figure. Indeed, 
the master equation (A2 is valid for non-monotonic con- 
tinuous functions r(t) even if there are discontinuities in 
its time derivative, provided dr/dt < 1, cf. lower panel 
of Fig. 11. 

In summary, as long as r(t) is a continuous function of 
time and dr/dt < 1 we find that Eq. (Al) holds, which in 
turn implies the master equation (A2). This is of course 
in general a complicated equation for P2(t), nonlinear 
and nonlocal in time. 

Next let us suppose that r(t) itself has discontinuities. 
Two generic situations are possible, as shown in Fig. 12. 
First consider a discontinuous abrupt decrease of r at t = 
t d from a high value of r_ to the left of the discontinuity 
to a low value r+ to the right, as seen in the upper panels. 
For later reference we call this Case 1. The behavior 
away from to on either side may be as described above, 
in which case it requires no extra analysis, or increases 
with time and violates the condition dr/dt < 1, a case 
discussed subsequently. In any case, we must be careful 
as we approach to ■ The sketch in the upper left panel of 
Fig. 12 shows that approaching to from below gives 

lim G(M') = -*d +r_). (A3) 

t^t D 

Approaching the discontinuity from above (right upper 
panel) leads to 

lim 6(£,t') = 6(t' -t D +r + ). (A4) 

Thus we see that inherits the discontinuity of r(t), but 
this is automatically "built in" in Eq. (Al), which is still 




t - x(tj t t - z(l) t 
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t,-x(t,) t, t 2 -z(t 2 ) h 



FIG. 11. Left upper panel: Schematic of the construct lead- 
ing to the generalized master equation Eq. (A2) when r(P2(t)) 
is analytic and decreasing with time. Right upper panel: The 
corresponding schematic when r(t) is analytic and increasing 
with time. The master equation (A2) again results provided 
dr/dt < 1. Lower panel: again Eq. (A2) applies for contin- 
uous r(t) even if there are discontinuities in its derivative, 
provided dr/dt < 1. 



valid. Hence the generalized master equation is again 
given by Eq. (A2). 

Now suppose the discontinuity occurs in the opposite 
direction, from a low value r_ to a high value t + as we 
move from left to right. This is shown in the lower pan- 
els of Fig. 12. We call this Case 2. Again, sufficiently 
far from the discontinuity either the previous discussion 
surrounding Fig. 11 continues to hold, or the subsequent 
discussion does. As we approach to from the left all con- 
tinues as before up to the value t' = to — t_, 

lim. Q(t,t') = 9(t' -t D +r_). (A5) 
t~>t D 

This is clarified in the left lower panel of Fig. 12. How- 
ever, when t' lies in the interval [to — T+, to — t_] then 
for t" in the range [i, io] lying between t' and io wc 
see from the figure that the inequality is not satisfied, 
i.e. that t' - t' > r(t"). Here t' - i = r(t). However, 
when i" continues to move up within the range [t',to], 
the inequality is once again satisfied, t' —t < r(t"). The 
change from "not satisfied" to "satisfied" occurs at the 
value t' = t o — t_ , so that 

lim e(t,t') = 6{t' - t D +t_). (A6) 

*->•*+ 

Comparing the last two equations, we see that is now in 
fact continuous in spite of the discontinuity in r(t), that 
is, now does not inherit the discontuity as it did when 
r(t) abruptly decreases. This asymmetry in behavior is 
of course a consequence of the direction of the arrow of 
time. 

Furthermore, additional "different" behavior is ob- 
served as we continue further upward in time. Suppose 
we now look at values of t in the next range, [to,t*], 
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FIG. 12. Cases that can arise when r(t) has discontinuities. 
Upper panels: high to low discontinuities. Lower panels: low 
to high discontinuities. 
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